rm(list=ls())
set.seed(08536)
if(!is.null(sessionInfo()$otherPkgs)){
  invisible(lapply(paste0('package:', names(sessionInfo()$otherPkgs)), detach, character.only=TRUE, unload=TRUE))
}
if (!require(pacman)){
  install.packages(pacman)
  library(pacman)
}
p_load(stargazer,
       survey,
       dplyr,
       questionr,
       srvyr,
       jtools,
       interactions,
       ggplot2,
       coefplot,
       corrplot,
       sjPlot,
       interplot,
       ggeffects,
       viridis,
       mice,
       labelled,
       mitools,
       dotwhisker,
       texreg,
       gtsummary,
       randomForest,
       caret,
       poliscidata)
extract.mi<-function(model){
  sink("file")
  s<-summary(model) 
  names<-names(model$coefficients) 
  co<-s$results
  se<-s$se
  pval<-2*pnorm(abs((co/se)),lower.tail = F)
  nimp <- model$nimp
  # rs<-s$r.squared 
  # adj<-s$adj.r.squared 
  # n<-nobs(model) 
  gof<-c(nimp) 
  gof.names<-c("Num.\\imputations") 
  tr<-createTexreg(
    coef.names=names, 
    coef=co, 
    se=se, 
    pvalues=pval, 
    gof.names=gof.names, 
    gof=gof) 
  sink()
  return(tr)
}
setMethod("extract",signature=className("MIresult","mitools"), definition=extract.mi)
dat <- readRDS("../Data/clean_data.RDS") %>% remove_labels()
dat$hist[is.na(dat$hist)] <- 99
#Not run
# factors <- dat %>% dplyr::select_if(is.factor)
# dat$territ_loss_concern_scaled <- as.vector(dat$territ_loss_concern_scaled)
# datimpute <- dat %>% dplyr::select_if(is.numeric) %>% mice(.,method="cart")
# datimpute <- as.mids(complete(datimpute,"long",include=T) %>% bind_cols(bind_rows(replicate(6,factors,simplify = F)))
# saveRDS(datimpute,"data_imputed.RDS")
svyimpute <- imputationList(complete(readRDS("../Data/data_imputed.RDS"),"all"))
dat4 <- readRDS("../Data/clean_data_noV4.RDS") %>% mutate(territ_loss_concern_scaled = scales::rescale(territ_loss_concern))

meplot <- function(model,var1,var2,ci=.95
                   ){
  alpha <- 1-ci
  z <- qnorm(1-alpha/2)
  beta.hat <- coef(model)
  cov <- vcov(model)
  z0 <- seq(0,100,length.out=101)
  dy.dx <- beta.hat[var1] + 
    beta.hat[(length(beta.hat)-1)]*z0 + 
    beta.hat[(length(beta.hat))]*z0^2
  se.dy.dx <- sqrt(cov[var1,var1] + 
                     z0^2*cov[(nrow(cov)-1),(ncol(cov)-1)] + 
                     z0^4*cov[nrow(cov),ncol(cov)] + 
                     2*z0*cov[var1,(ncol(cov)-1)] +
                     2*(z0^2)*cov[var1,ncol(cov)] +
                     2*(z0^3)*cov[(nrow(cov)-1),ncol(cov)]
                   )
  upr <- dy.dx + z*se.dy.dx
  lwr <- dy.dx - z*se.dy.dx
  return(data.frame(x=z0,est=dy.dx,ci.upper=upr,ci.lower=lwr))
}
makeform <- function(country,outcome,controls=T,interact = F, obs=F,diaspora=F,fe=T,scale=F,quadratic=F){
  #Edit to change baseline controls
  treat <- case_when(diaspora ~ "D",
                     obs&(!scale) ~ "territ_loss_concern",
                     obs&scale ~ "territ_loss_concern_scaled",
                     TRUE ~ "territorial")
  if(interact){
    if(!quadratic){
      basespec <-c(paste(ifelse(scale,"territ_loss_concern_scaled","territ_loss_concern"),treat,sep="*"), 
                   "age", 
                   "educ_univ_complete", 
                   "educ_sec_univprep", 
                   "educ_sec_techvoc", 
                   "male", 
                   "city",
                   "town",
                   "covid_civlib_restrict_ind",
                   "econ_intervention_ind",
                   "pro_global_ind",
                   "democ_ind",
                   "natid_ind",
                   "neigh_lgbt",
                   "relig_import01",
                   "pol_int01")
    } else{
      basespec <-c(paste(ifelse(scale,"territ_loss_concern_scaled","territ_loss_concern"),treat,sep="*"),
                   paste(ifelse(scale,"I(territ_loss_concern_scaled^2)","I(territ_loss_concern^2)"),treat,sep="*"),
                   "age", 
                   "educ_univ_complete", 
                   "educ_sec_univprep", 
                   "educ_sec_techvoc", 
                   "male", 
                   "city",
                   "town",
                   "covid_civlib_restrict_ind",
                   "econ_intervention_ind",
                   "pro_global_ind",
                   "democ_ind",
                   "natid_ind",
                   "neigh_lgbt",
                   "relig_import01",
                   "pol_int01")
    }
    
  } else{
    basespec <-c(treat,
                 "age", 
                 "educ_univ_complete", 
                 "educ_sec_univprep", 
                 "educ_sec_techvoc", 
                 "male", 
                 "city",
                 "town",
                 "covid_civlib_restrict_ind",
                 "econ_intervention_ind",
                 "pro_global_ind",
                 "democ_ind",
                 "natid_ind",
                 "neigh_lgbt",
                 "relig_import01",
                 "pol_int01")
  }
  if(!controls){
    if(interact){
      if(!quadratic){
        basespec <- paste(ifelse(scale,"territ_loss_concern_scaled","territ_loss_concern"),treat,sep="*")
      } else{
        basespec <-c(paste(ifelse(scale,"territ_loss_concern_scaled","territ_loss_concern"),treat,sep="*"),
                     paste(ifelse(scale,"I(territ_loss_concern_scaled^2)","I(territ_loss_concern^2)"),treat,sep="*")
                     )
      }
    } else{
      basespec <- treat
    }
  }
  #Add country-appropriate region fixed effects
  
    geo <- case_when(country == "Romania" ~ "County_Rom",
                     country == "Hungary" ~ "County_Hung",
                     country == "Turkey" ~ "Province_Turk",
                     country == "Germany" ~ "Bundesland_Germ",
                     country == "All" ~ "Country"
    )
  
  
  
  if(outcome == "territ_loss_concern"){
      basespec <-c("age", 
                   "educ_univ_complete", 
                   "educ_sec_univprep", 
                   "educ_sec_techvoc",
                   "male", 
                   "national_language",
                   "covid_civlib_restrict_ind",
                   "econ_intervention_ind",
                   "pro_global_ind",
                   "democ_ind",
                   "neigh_lgbt",
                   "relig_import01",
                   "pol_int01",
                   "cultid",
                   "pplid",
                   "statid",
                   "town"
                   )
      if(country == "Hungary"){
        polit <- c("like_jobbik", 
                   "like_mszp",
                   "like_fidesz")
      } else if(country == "Romania"){
        polit <- c("like_prm",
                   "like_psd",
                   "like_pmp",
                   "like_udmr",
                   "like_usr",
                   "like_pnl")
      }else if(country == "Turkey"){
        polit <- c("like_akp", # Turkey
                   "like_mhp",
                   "like_hdp",
                   "like_chp",
                   "like_iyi")
      } else if(country == "Germany"){
        polit <- c("like_cducsu", 
                   "like_spd",
                   "like_afd")
      } else{
        polit <- c("nationalist_like_all_onep")
      }
      basespec <- c(basespec,polit)
  }

  f <- as.formula(
    paste(outcome, 
          ifelse(fe,paste(c(basespec,geo), collapse = " + "),paste(basespec, collapse = " + ")), 
          sep = " ~ "))
  return(f)
  
}
getmodels <- function(outcomes,countries,binary=F,controls=T,interact = F, obs=F,diaspora=F,fe=T,scale=F,quadratic=F){
  models <- list()
  if(binary){
    for(i in seq_along(outcomes)){
      if(countries[i] == "Hungary"){
        models[[i]] <- hung %>% svyglm(makeform(countries[i],outcomes[i],controls,interact, obs,diaspora,fe,scale,quadratic),.,family = quasibinomial())
      } else if(countries[i] == "Romania"){
        models[[i]] <- roma %>% svyglm(makeform(countries[i],outcomes[i],controls,interact, obs,diaspora,fe,scale,quadratic),.,family = quasibinomial())
      }else if(countries[i] == "Turkey"){
        models[[i]] <- turk %>% svyglm(makeform(countries[i],outcomes[i],controls,interact, obs,diaspora,fe,scale,quadratic),.,family = quasibinomial())
      } else if(countries[i] == "Germany"){
        models[[i]] <- germ %>% svyglm(makeform(countries[i],outcomes[i],controls,interact, obs,diaspora,fe,scale,quadratic),.,family = quasibinomial())
      } else{
        models[[i]] <- svydat %>% svyglm(makeform(countries[i],outcomes[i],controls,interact, obs,diaspora,fe,scale,quadratic),.,family = quasibinomial())
      }
    }
  } else if (outcomes[[1]] == "territ_loss_concern"){
    for(i in seq_along(outcomes)){
      if(countries[i] == "Hungary"){
        models[[i]] <- hungall %>% svyglm(makeform(countries[i],outcomes[i],controls,interact, obs,diaspora,fe,scale,quadratic),.)
      } else if(countries[i] == "Romania"){
        models[[i]] <- romaall %>% svyglm(makeform(countries[i],outcomes[i],controls,interact, obs,diaspora,fe,scale,quadratic),.)
      }else if(countries[i] == "Turkey"){
        models[[i]] <- turkall %>% svyglm(makeform(countries[i],outcomes[i],controls,interact, obs,diaspora,fe,scale,quadratic),.)
      } else if(countries[i] == "Germany"){
        models[[i]] <- germall %>% svyglm(makeform(countries[i],outcomes[i],controls,interact, obs,diaspora,fe,scale,quadratic),.)
      } else{
        models[[i]] <- svydatall %>% svyglm(makeform(countries[i],outcomes[i],controls,interact, obs,diaspora,fe,scale,quadratic),.)
      }
    }
  } else{
    for(i in seq_along(outcomes)){
      if(countries[i] == "Hungary"){
        models[[i]] <- hung %>% svyglm(makeform(countries[i],outcomes[i],controls,interact, obs,diaspora,fe,scale,quadratic),.)
      } else if(countries[i] == "Romania"){
        models[[i]] <- roma %>% svyglm(makeform(countries[i],outcomes[i],controls,interact, obs,diaspora,fe,scale,quadratic),.)
      }else if(countries[i] == "Turkey"){
        models[[i]] <- turk %>% svyglm(makeform(countries[i],outcomes[i],controls,interact, obs,diaspora,fe,scale,quadratic),.)
      } else if(countries[i] == "Germany"){
        models[[i]] <- germ %>% svyglm(makeform(countries[i],outcomes[i],controls,interact, obs,diaspora,fe,scale,quadratic),.)
      } else{
        models[[i]] <- svydat %>% svyglm(makeform(countries[i],outcomes[i],controls,interact, obs,diaspora,fe,scale,quadratic),.)
      }
    }
  }
  
  return(models)
}
svydat <-dat %>% filter(national_language==1) %>% svydesign(~1, strata = ~Country, weights=~wt,data=.)
svydatimpute <- subset(svydesign(~1, strata = ~Country, weights=~wt,data=svyimpute),national_language==1)
svydatimputehuntur <- subset(svydatimpute,hun==1|tur==1)
svydat4impute <-  subset(svydatimpute,Version!=4)
svydat4imputehuntur <- subset(svydatimputehuntur,Version!=4)
svydat4 <-  subset(svydat,Version!=4)
svydat4votedpp <-  subset(svydat4,populist_voted==1)
svydat4novotepp <-  subset(svydat4,populist_voted==0)
svydatvotedpp <-  subset(svydat,populist_voted==1)
svydatnovotepp <-  subset(svydat,populist_voted==0)
svydatall <-dat %>% svydesign(~1, strata = ~Country, weights=~wt,data=.)

# Table F.1
t <- tbl_svysummary(svydat, include= c(
  "nationalist_like_all_max",
  "populist_voted",
  "hist_greatestextent",
  "territ_loss_concern",
  "age", 
  "educ_univ_complete", 
  "educ_sec_univprep", 
  "educ_sec_techvoc", 
  "male", 
  "city",
  "town",
  "covid_civlib_restrict_ind",
  "econ_intervention_ind",
  "pro_global_ind",
  "democ_ind",
  "natid_ind",
  "neigh_lgbt",
  "relig_import01",
  "pol_int01"),
  statistic=list(all_continuous() ~ "& {mean} & {sd} & {min} & {max} \\\\", all_categorical() ~ "0.{p}"), missing = "no",
  type = list(everything() ~ "continuous"))
t

hung <- dat %>% filter(hun==1,!is.na(County_Hung),national_language==1) %>% svydesign(~1,weights=~wt,data = .)
roma <- dat %>% filter(rom==1,!is.na(County_Rom),national_language==1) %>% svydesign(~1,weights=~wt,data = .)
turk <- dat %>% filter(tur==1,!is.na(Province_Turk),national_language==1) %>% svydesign(~1,weights=~wt,data = .)
germ <- dat %>% filter(ger==1,!is.na(Bundesland_Germ),national_language==1) %>% svydesign(~1,weights=~wt,data = .)
hungimpute <- subset(svydatimpute,hun==1)
romaimpute <- subset(svydatimpute,rom==1)
turkimpute <- subset(svydatimpute,tur==1)
germimpute <- subset(svydatimpute,ger==1)

hung4 <- dat4 %>% filter(hun==1,!is.na(County_Hung),national_language==1) %>% svydesign(~1,weights=~wt,data = .)
roma4 <- dat4 %>% filter(rom==1,!is.na(County_Rom),national_language==1) %>% svydesign(~1,weights=~wt,data = .)
turk4 <- dat4 %>% filter(tur==1,!is.na(Province_Turk),national_language==1) %>% svydesign(~1,weights=~wt,data = .)
germ4 <- dat4 %>% filter(ger==1,!is.na(Bundesland_Germ),national_language==1) %>% svydesign(~1,weights=~wt,data = .)

hung4impute <- subset(svydat4impute,hun==1)
roma4impute <- subset(svydat4impute,rom==1)
turk4impute <- subset(svydat4impute,tur==1)
germ4impute <- subset(svydat4impute,ger==1)

hungall <- dat %>% filter(hun==1,!is.na(County_Hung)) %>% svydesign(~1,weights=~wt,data = .)
romaall <- dat %>% filter(rom==1,!is.na(County_Rom)) %>% svydesign(~1,weights=~wt,data = .)
turkall <- dat %>% filter(tur==1,!is.na(Province_Turk)) %>% svydesign(~1,weights=~wt,data = .)
germall <- dat %>% filter(ger==1,!is.na(Bundesland_Germ)) %>% svydesign(~1,weights=~wt,data = .)

#########################################################################
###############Main Models -- nationalists ##############################
#########################################################################

#Table A.1
models <- list(
  with(svydatimpute,svyglm(makeform("All","nationalist_like_all_max",
                  obs=T,
                  controls=F,
                  fe=F,
                  scale = T))),
  with(svydatimpute,svyglm(makeform("All","nationalist_like_all_max",
                  obs=T,
                  controls=F,
                  fe=T,
                  scale = T))),
  with(svydat4impute,svyglm(makeform("All","nationalist_like_all_max",
                                     obs=F,
                                     controls=F,
                                     fe=F))),
  with(svydat4impute,svyglm(makeform("All","nationalist_like_all_max",
                                     obs=F,
                                     controls=T,
                                     fe=T))),
  with(svydat4impute,svyglm(makeform("All","nationalist_like_all_max",
                                     interact=T,
                                     controls=F,
                                     fe=F,
                                     scale = T))),
  with(svydat4impute,svyglm(makeform("All","nationalist_like_all_max",
                                     interact=T,
                                     controls=T,
                                     fe=T,
                                     scale = T)))
)

out <- lapply(models,MIcombine)
rsquared <- round(unlist(lapply(models,function(x){mean(sapply(x,fit.svyglm)[1,])})),2)
paste0(c("R$^2$",rsquared),collapse = " & ")

texreg(out,
       omit.coef=c('Country'),
       custom.coef.names = c("Constant",
                             "Territorial Loss Concern",
                             "Territorial Prime",
                             "Age",
                             "Completed University",
                             "Completed Secondary",
                             "Technical/Vocational Education",
                             "Male",
                             "Urban",
                             "Small Town",
                             "COVID Restriction Support Index",
                             "Interventionism Index",
                             "Globalization Support Index",
                             "Support for Democracy Index",
                             "National Identity Index",
                             "Accept LGBT Neighbors",
                             "Religious Importance",
                             "Political Interest",
                             "Prime : Concern"),
       stars=c(0.01,0.05,0.1)
)

#########################################################################
###############die Linke Models -- Table A.2#############################
#########################################################################


models <- list(
  with(germimpute,svyglm(makeform("Germany","populist_only_ger",
                                            obs=T,
                                            controls=F,
                                            fe=F,
                                            scale=T))),
  with(germimpute,svyglm(makeform("Germany","populist_only_ger",
                                            obs=T,
                                            controls=T,
                                            fe=T,
                                            scale=T))),
  with(germ4impute,svyglm(makeform("Germany","populist_only_ger",
                                             controls=F,
                                             fe=F,
                                             scale=T))),
  with(germ4impute,svyglm(makeform("Germany","populist_only_ger",
                                             controls=T,
                                             fe=T,
                                             scale=T))),
  with(germ4impute,svyglm(makeform("Germany","populist_only_ger",interact = T,
                                             quadratic=F,
                                             controls=F,
                                             fe = F,
                                             scale=T))),
  with(germ4impute,svyglm(makeform("Germany","populist_only_ger",interact = T,
                                             quadratic=F,
                                             controls=T,
                                             fe = T,
                                             scale=T)))
)

out <- lapply(models,MIcombine)
rsquared <- round(unlist(lapply(models,function(x){mean(sapply(x,fit.svyglm)[1,])})),2)
paste0(c("R$^2$",rsquared),collapse = " & ")
texreg(out,
       omit.coef=c('Bundesland'),
       custom.coef.names = c("Constant",
                             "Territorial Loss Concern",
                             "Age",
                             "Completed University",
                             "Completed Secondary",
                             "Technical/Vocational Education",
                             "Male",
                             "Urban",
                             "Small Town",
                             "COVID Restriction Support Index",
                             "Interventionism Index",
                             "Globalization Support Index",
                             "Support for Democracy Index",
                             "National Identity Index",
                             "Accept LGBT Neighbors",
                             "Religious Importance",
                             "Political Interest",
                             "Territorial Prime",
                             "Prime : Concern"),
       stars=c(0.01,0.05,0.1)
)
#########################################################################
###############Hungary Lost Unity Only Models -- Table B.1###############
#########################################################################


models <- list(
  with(hungimpute,svyglm(makeform("Hungary","populist_like_all_max",
                                            obs=T,
                                            controls=F,
                                            fe=F,
                                            scale=T))),
  with(hungimpute,svyglm(makeform("Hungary","populist_like_all_max",
                                            obs=T,
                                            controls=F,
                                            fe=T,
                                            scale=T))),
  with(hungimpute,svyglm(makeform("Hungary","populist_like_all_max",
                                            obs=T,
                                            controls=T,
                                            fe=T,
                                            scale=T)))
)
out <- lapply(models,MIcombine)
rsquared <- round(unlist(lapply(models,function(x){mean(sapply(x,fit.svyglm)[1,])})),2)
paste0(c("R$^2$",rsquared),collapse = " & ")
texreg(out,
       omit.coef=c('County|Bundesland|Province'),
       custom.coef.names = c("Constant",
                             "Territorial Loss Concern",
                             "Age",
                             "Completed University",
                             "Completed Secondary",
                             "Technical/Vocational Education",
                             "Male",
                             "Urban",
                             "Small Town",
                             "COVID Restriction Support Index",
                             "Interventionism Index",
                             "Globalization Support Index",
                             "Support for Democracy Index",
                             "National Identity Index",
                             "Accept LGBT Neighbors",
                             "Religious Importance",
                             "Political Interest"),
       stars=c(0.01,0.05,0.1)
)

#########################################################################
###############Romania Models -- Table B.2  #############################
#########################################################################
romaterritory <- subset(romaimpute,territoryonly==1)
romalost <- subset(romaimpute,lostunity==1)

models <- list(
  with(romaterritory,svyglm(makeform("Romania","populist_like_all_max",
                                               obs=T,
                                               controls=F,
                                               fe=F,
                                               scale=T))),
  with(romaterritory,svyglm(makeform("Romania","populist_like_all_max",
                                               obs=T,
                                               controls=F,
                                               fe=T,
                                               scale=T))),
  with(romaterritory,svyglm(makeform("Romania","populist_like_all_max",
                                               obs=T,
                                               controls=T,
                                               fe=T,
                                               scale=T))),
  with(romalost,svyglm(makeform("Romania","populist_like_all_max",
                                     obs=T,
                                     controls=F,
                                     fe=F,
                                     scale=T))),
  with(romalost,svyglm(makeform("Romania","populist_like_all_max",
                                     obs=T,
                                     controls=F,
                                     fe=T,
                                     scale=T))),
  with(romalost,svyglm(makeform("Romania","populist_like_all_max",
                                     obs=T,
                                     controls=T,
                                     fe=T,
                                     scale=T)))
)
out <- lapply(models,MIcombine)
rsquared <- round(unlist(lapply(models,function(x){mean(sapply(x,fit.svyglm)[1,])})),2)
paste0(c("R$^2$",rsquared),collapse = " & ")
texreg(out,
       omit.coef=c('County|Bundesland|Province'),
       custom.coef.names = c("Constant",
                             "Territorial Loss Concern",
                             "Age",
                             "Completed University",
                             "Completed Secondary",
                             "Technical/Vocational Education",
                             "Male",
                             "Urban",
                             "Small Town",
                             "COVID Restriction Support Index",
                             "Interventionism Index",
                             "Globalization Support Index",
                             "Support for Democracy Index",
                             "National Identity Index",
                             "Accept LGBT Neighbors",
                             "Religious Importance",
                             "Political Interest"),
       stars=c(0.01,0.05,0.1)
)



#########################################################################
###############Germany Territory Only Models -- Table B.3################
#########################################################################
germterritory <- subset(germimpute,territoryonly==1)
germlost <- subset(germimpute,lostunity==1)
models <- list(
  with(germterritory,svyglm(makeform("Germany","populist_like_all_max",
                                               obs=T,
                                               controls=F,
                                               fe=F,
                                               scale=T))),
  with(germterritory,svyglm(makeform("Germany","populist_like_all_max",
                                               obs=T,
                                               controls=F,
                                               fe=T,
                                               scale=T))),
  with(germterritory,svyglm(makeform("Germany","populist_like_all_max",
                                               obs=T,
                                               controls=T,
                                               fe=T,
                                               scale=T))),
  with(germlost,svyglm(makeform("Germany","populist_like_all_max",
                                     obs=T,
                                     controls=F,
                                     fe=F,
                                     scale=T))),
  with(germlost,svyglm(makeform("Germany","populist_like_all_max",
                                     obs=T,
                                     controls=F,
                                     fe=T,
                                     scale=T))),
  with(germlost,svyglm(makeform("Germany","populist_like_all_max",
                                     obs=T,
                                     controls=T,
                                     fe=T,
                                     scale=T)))
)
out <- lapply(models,MIcombine)
rsquared <- round(unlist(lapply(models,function(x){mean(sapply(x,fit.svyglm)[1,])})),2)
paste0(c("R$^2$",rsquared),collapse = " & ")
texreg(out,
       omit.coef=c('County|Bundesland|Province'),
       custom.coef.names = c("Constant",
                             "Territorial Loss Concern",
                             "Age",
                             "Completed University",
                             "Completed Secondary",
                             "Technical/Vocational Education",
                             "Male",
                             "Urban",
                             "Small Town",
                             "COVID Restriction Support Index",
                             "Interventionism Index",
                             "Globalization Support Index",
                             "Support for Democracy Index",
                             "National Identity Index",
                             "Accept LGBT Neighbors",
                             "Religious Importance",
                             "Political Interest"),
       stars=c(0.01,0.05,0.1)
)

#########################################################################
###############Turkey Territory Only Models -- Table B.4#################
#########################################################################
turkterritory <- subset(turkimpute,territoryonly==1)
turklost <- subset(turkimpute,lostunity==1)
models <- list(
  with(turkterritory,svyglm(makeform("Turkey","populist_like_all_max",
                                               obs=T,
                                               controls=F,
                                               fe=F,
                                               scale=T))),
  with(turkterritory,svyglm(makeform("Turkey","populist_like_all_max",
                                               obs=T,
                                               controls=F,
                                               fe=T,
                                               scale=T))),
  with(turkterritory,svyglm(makeform("Turkey","populist_like_all_max",
                                               obs=T,
                                               controls=T,
                                               fe=T,
                                               scale=T))),
  with(turklost,svyglm(makeform("Turkey","populist_like_all_max",
                                     obs=T,
                                     controls=F,
                                     fe=F,
                                     scale=T))),
  with(turklost,svyglm(makeform("Turkey","populist_like_all_max",
                                     obs=T,
                                     controls=F,
                                     fe=T,
                                     scale=T))),
  with(turklost,svyglm(makeform("Turkey","populist_like_all_max",
                                     obs=T,
                                     controls=T,
                                     fe=T,
                                     scale=T)))
)
out <- lapply(models,MIcombine)
rsquared <- round(unlist(lapply(models,function(x){mean(sapply(x,fit.svyglm)[1,])})),2)
paste0(c("R$^2$",rsquared),collapse = " & ")
texreg(out,
       omit.coef=c('County|Bundesland|Province'),
       custom.coef.names = c("Constant",
                             "Territorial Loss Concern",
                             "Age",
                             "Completed University",
                             "Completed Secondary",
                             "Technical/Vocational Education",
                             "Male",
                             "Urban",
                             "Small Town",
                             "COVID Restriction Support Index",
                             "Interventionism Index",
                             "Globalization Support Index",
                             "Support for Democracy Index",
                             "National Identity Index",
                             "Accept LGBT Neighbors",
                             "Religious Importance",
                             "Political Interest"),
       stars=c(0.01,0.05,0.1)
)

#########################################################################
###############Turkey Territory Only Models -- Kemalists, Table B.5######
#########################################################################

models <- list(
  with(turkterritory,svyglm(makeform("Turkey","nationalist_like_all_max",
                                               obs=T,
                                               controls=F,
                                               fe=F,
                                               scale=T))),
  with(turkterritory,svyglm(makeform("Turkey","nationalist_like_all_max",
                                               obs=T,
                                               controls=F,
                                               fe=T,
                                               scale=T))),
  with(turkterritory,svyglm(makeform("Turkey","nationalist_like_all_max",
                                               obs=T,
                                               controls=T,
                                               fe=T,
                                               scale=T))),
  with(turklost,svyglm(makeform("Turkey","nationalist_like_all_max",
                                     obs=T,
                                     controls=F,
                                     fe=F,
                                     scale=T))),
  with(turklost,svyglm(makeform("Turkey","nationalist_like_all_max",
                                     obs=T,
                                     controls=F,
                                     fe=T,
                                     scale=T))),
  with(turklost,svyglm(makeform("Turkey","nationalist_like_all_max",
                                     obs=T,
                                     controls=T,
                                     fe=T,
                                     scale=T)))
)
out <- lapply(models,MIcombine)
rsquared <- round(unlist(lapply(models,function(x){mean(sapply(x,fit.svyglm)[1,])})),2)
paste0(c("R$^2$",rsquared),collapse = " & ")
texreg(out,
       omit.coef=c('County|Bundesland|Province'),
       custom.coef.names = c("Constant",
                             "Territorial Loss Concern",
                             "Age",
                             "Completed University",
                             "Completed Secondary",
                             "Technical/Vocational Education",
                             "Male",
                             "Urban",
                             "Small Town",
                             "COVID Restriction Support Index",
                             "Interventionism Index",
                             "Globalization Support Index",
                             "Support for Democracy Index",
                             "National Identity Index",
                             "Accept LGBT Neighbors",
                             "Religious Importance",
                             "Political Interest"),
       stars=c(0.01,0.05,0.1)
)



